Quantum melting of a crystal of dipolar bosons 
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I. INTRODUCTION 

Quantum melting (the melting of a crystal induced 
by quantum fluctuations) can be very different from a 
classical melting (induced by temperature). In two di- 
mensional electronic systems, it was found that near the 
melting points, the system is neither in a liquid nor in 
a solid phase but shares properties with both phases.^ 
In bosonic systems, such an intermediate phase ("super- 
solid" between the solid and the superfluid phase) was 
proposed many years ago4££ The supersolid phase can 
be thought as a condensate of either static or dynamical 
defects^ Some supersolid behaviour (i.e. existence of a 
non zero superfluid fraction in a solid system) has indeed 
been recently observed^ in Helium 4, creating a large 
renew of interest on both theoretical and experimental 
sides. The situation in Helium 4 is however controversial 
as (i) uptodate quantum Monte-Carlo simulations do not 
find any superfluid fraction in the solid phased ' 10 : 11 while 
providing very precise predictions in the superfluid and 
the solid phase, (ii) the experiments seem to be sensi- 
tive to the way the solid is prepared^ (an annealing of 
the supersolid kills or enhances the superfluid fraction 
depending on the experiments) indicating that disorder 
probably plays a key role in this physics and (iii) one 
gets only access experimentally to macroscopic quantities 
which make it difficult to infer the physical mechanisms 
actually involved. 

An alternative bosonic system is now promising due 
to recent progresses in the field of ultracold gases. 
The production and observation of heteronuclear alkali 
molecules^ at ultralow temperatures has been reached 
with various dimers, e.g. KRb 14,15 , RbCs 16 or LiCs^, 
yet not in their true molecular ground state. Polar 



molecules confined in a two-dimensional trap with a per- 
pendicular electric field could be given an important 
dipole moment (on the order of 1 Debysi&) leading to 
strong repulsive dipole-dipole interaction on the plane. 
Although these interactions are not trully long range in 
2D, the potential range is larger than for typical van der 
Waals interaction and smaller densities are necessary to 
observe crystallization. Moreover the true repulsive char- 
acter of dipole-dipole interaction should strongly inhibit 
inelastic three-body recombination. Such systems could 
be very good candidates for the study of quantum melt- 
ing as they do not share the usual limitations found in 
Helium 4: there is a complete control of the external po- 
tential seen by molecules (no disorder) and many direct 
observations can been done (of the spatial distribution 
of the density for the solid, of the condensate fraction 
for the superfluid part). It is the purpose of this paper 
to study this potential use of cold dipolar gases in some 
detail. We note that two papers with a similar topic 
and some overlapping results appeared recentl y 19 : 20 and 
hence, whenever appropriate we will compare our find- 
ings with these works. Another possible experimental re- 
alization of the model studied in this paper are electron- 
hole heterostructure bilayers where the excitons also in- 
teract through dipole-dipole interaction^ 

After introducing the model (Bosons in two dimensions 
with a dipole-dipole repulsion), and briefly the numeri- 
cal technique (Green Function Monte-Carlo or GFMC) 
in section [Til we proceed with the description of the high 
density limit where the crystal is expected. In section Hill 
we calculate the phonon spectrum of the crystal, from 
which we extract the high density expansion of the energy 
of the system. The nature (first order) of the transition is 
then discussed in section HV1 using the GFMC technique. 
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Particular attention is paid to the sensitivity of the re- 
sults to the trial wave-function used in the calculation. 
The quality of the wave-function is usually measured in 
term of the distance in energy to the true ground state. 
We find that a better measure can be done by calculat- 
ing the overlap of the trial wave-function with the true 
ground-state^,. In section [V] we discuss the possibility 
of the existence of a supersolid in this system. In the last 
section I VII we estimate the regime of parameters that 
should be accessible experimentally. The appendix con- 
tains the Ewald summation technique that is applied to 
the dipole-dipole interaction to reduce finite size effects. 

II. THE MODEL: N BOSONS IN 2 DIMENSIONS 
WITH DIPOLE-DIPOLE INTERACTION 

Our system is made of N bosons of mass Mq confined in 
two dimensions which interact through a dipole-dipole in- 
teraction of strength Cdd, so that the Hamiltonian reads 

2—1 1<J ' J 1 

with dipoles polarized perpendicular to the plane. Mea- 
suring lengths in unit of the average distance I between 
bosons (defined as 7r^ 2 = 1/n where n is the 2D den- 
sity), and the energies in unit of Eg = Cdd/{2£ 3 ), the 
rescaled Hamiltonian depends on a unique dimensionless 
parameter r s and reads, 

1 N 1 

2 — 1 i<J ' J ' 

The parameter r s (named in analogy with electronic sys- 
tems) roughly measures the ratio between interaction 
energy ~ Cdd/(? and kinetic energy ~ h 2 /(Mo£ 2 ) and 
reads, 

r s = t B /i, (3) 

where t* B = MoCdd/Ji 2 - Contrary to electronic systems, 
r s increases with density, so that Bose-Einstein conden- 
sation is obtained at low density while at high density 
the dipole-dipole interaction dominates and the system 
crystallizes into a triangular lattice (see Sec. IIIip . 

A. The Green Function Monte-Carlo technique 

The ground state properties of the Hamiltonian Eq. ([2]) 
can be conveniently studied using quantum Monte-Carlo 
techniques which for a bosonic system (in contrast to 
fcrmionic systems where one has to deal with the so called 
sign problem) provides access to the exact ground state 
|\&o) °f the system. The idea behind this method is to 
project an initial (variational) Guiding Wave Function 
I^g) (GWF) on the exact ground state \^>o) by applying 



the projection operator e~@ H . For large projecting time 
/?, the excited states are exponentially filtered out from 
the initial trial wave- function, and one gets, 

|* ) oc lim e~ 0H \^ a ) (4) 

GFMC is a scheme to implement Eq. Q in a stochastic 
way. In practice, |\&(?) n ot only enters as the starting 
point but also through the way the Hilbert space is sam- 
pled when the projection operator is applied. A good 
choice of the trial wave-function is therefore important 
to get a good signal to noise ratio. In principle however 
(i.e. for infinite computing time), the results should be 
independent of Y&g) ■ 

Although this technique can be used for continuum 
systems, our particular implementation (which has been 
described in Ref. y) relies on the discretization of the 
Hamiltonian on a grid (which in returns allows us to work 
without discretizing the imaginary time f3) . For some cal- 
culations, we have also used the recently introduced Rep- 
tation Quantum Monte Carlo algorithm (RQMC)22i2£. 
RQMC is formally equivalent to GFMC and is its path 
integral counter part. For the regimes discussed in this 
paper, we found however that RQMC gave no real ad- 
vantage over GFMC. Unless explicitly stated we use only 
pure estimators (forward walking technique^) for the 
quantum averages. 

B. Discretized Hamiltonian 

The bosonic system is put on a rectangular grid of 
L x x L y sites with nearest neighbor hopping. Finite size 
effects are considerably reduced by repeating periodically 
the simulation box on the two-dimensional plane. This 
amounts to use periodic boundary conditions for hop- 
ping and to replace the dipole interaction by an effective 
two-body interaction that includes all periodic images. 
The resulting summation is efficiently performed using a 
straightforward extension of the Ewald summation tech- 
nique. This is discussed in more details in appendix [X] 
The discretized Hamiltonian on the grid reads 

H = -t^2 etc?, + j^2v(f- f')npn P , + {At + A) AT, 

{r,r') r^r' 

(5) 

where the operator ct (cp) creates (destroys) a boson on 
point r with the standard commutation relation rules, 
the sum r') ^ s done on the nearest neighbor points 
on the grid. The form of the periodized two-body inter- 
action V(r) and of the constant A (which accounts for the 
interaction of a boson with its own images) are given in 
appendix[X] Noting v = N/(L x L y ) the filling factor, the 
continuous model Eq.® is reproduced by choosing the 
hopping amplitude to be t = l/{r s T:v) and the effective 
interaction strength U = 2/(nv) 3 / 2 in the limit v —* . 
In practice we took v — 1/56 (some small influence of the 
presence of the grid) down to v — 1/780 (no detectable 
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influence of the grid). The grid dimensions are also cho- 
sen in order not to induce distortion of the triangular 
lattice of the crystal. 

C. The guiding wave- functions 

All our Monte Carlo simulations start from a trial 
wavefunction which is an initial guess for the exact 
ground state wavefunction. The trial GWF is also used 
as a guide to converge to the ground state during the 
GFMC projection. For practical efficiency, i.e. to reduce 
noise in numerical results, it is important that the GWF 
is as close as possible from the unknown exact solution. 
We use here a Bijl-Jastrow form: 

N 

*GWF(ri, ■ ,T N ) = Y[ M*i) II M r 3 ~ rfe D ( 6 ) 

i—1 j<k 

which includes at the variational level one-body and two- 
body correlations. Two-body correlations are essential 
in the homogeneous (BEC) limit where it dominates the 
physics and determines the shape of the wavefunction. 
On the contrary, one-body terms are sufficient to de- 
scribe the crystal limit where atoms are almost frozen. It 
is worth noticing that the two-body scattering problem 
can be solved exactly at zero energy, leading </>2 _b (r) = 
AKo(2y/r s /r), where A is an arbitrary constant and Kq 
the modified Besscl function of the second kind. Its s hort 
distance behaviour is 2 ~ b (r) cx (r /r s ) 1 / A e~' 2 ^ r ^ r . In 
practice, we chose 

■fc(r) = exp {- 2 J±e-r'A, (7) 

with a single variational parameter A which sets the 
length over which two-body correlations are suppressed. 
In particular, it reproduces the leading short distance be- 
havior for the two-body scattering problem. More com- 
plicated variational forms have been tested (with more 
variational parameters) with the same efficiency. Typi- 
cal values of A ~ 1 have been found in the vicinity of the 
quantum melting point r s ~ 27. 

Noting Ay (Ax) the distance between sites along 
the y (x) axis (Ay /Ax = y/3/2), we define the vec- 
tor qi = (0, 2%/ Ay) in the reciprocal lattice. Vectors 
q 2 /3 = (±2n/Ax, — n/Ay) are obtained by rotating qi 
by an angle of 2tt/3 and —2ir/3. We choose for the one- 
body trial function 

M r ) = II (1 + " cos( qi •!■)), (8) 

i=l,2,3 

where maxima reproduce the triangular lattice expected 
for the crystal. This function is well-suited to describe 
quantum melting since it interpolates between a flat liq- 
uid (BEC)-type pattern for a = and a triangular crystal 
form for a =/= 0. Moreover it conserves the symmetriza- 
tion property imposed by bosonic statistics in contrast 



with the Gaussian form often considered in the litera- 
ture. 



III. PROPERTIES OF THE CRYSTAL PHASE 

The infinite r s limit corresponds to a purely classi- 
cal system of dipoles with interaction energy only. We 
have compared energies of various lattice configurations 
and found that the ground state is given by the tri- 
angular Bravais lattice of energy per particle E$ — 
1.59702.... For comparison, the square lattice energy 
is E — 1.62232 . . .. Ground state properties of the crys- 
tal phase can be computed perturbatively for large r s 
following Refi2£. For large but finite r s , kinetic energy 
plays a role as atoms are able to move around the perfect 
triangular lattice configuration. We treat kinetic terms 
perturbatively which amounts to an harmonic approxi- 
mation for the potential felt by the atoms. The emerging 
collective modes are then a sum of harmonic oscillators 
over the Brillouin zone. They describe phonons and dom- 
inate dynamical properties of the ground state. Moreover 
first order correction in \jyfrl to the ground state en- 
ergy derives from zero-point energies of these vibrational 
modes. 

The Hamiltonian Eq. §2$ expands around the classical 
stable position r$ = Ri + where {Ri} denotes the tri- 
angular lattice configuration. Stopping at second order in 
atomic displacements, it is straightforward to diagonalize 
the Hamiltonian in Fourier space. Collective eigenmodes 
read Uj(i) = u q e^ q ' Rj _wt ^ where the vectors u q are so- 
lutions of the eigenvalue problem 

w 2 (q)u q = -*(q)flq. (9) 

r s 

Here $(q) is a tensor or 2 x 2 matrix defined by 

u (i°o) 

with Cartesian coordinates a,/3 — x,y. For each 
wavevector q in the first Brillouin zone, Eq. ([9]) has two 
solutions which leads to two branches of collective exci- 
tations. The symmetry of the triangular lattice under 
the point group Cq v allows to restrict values of q within 
the irreducible Brillouin zone represented in FigfTJ For 
illustration we show in Fig. [1] the corresponding phonon 
eigenmodes, solutions of Eq. @ , along the boundary of 
the irreducible Brillouin zone. We rescale Eq. ([9]) with 
the definition 

where b — J 2j\fi is the lattice spacing for the trian- 
gular lattice with density l/ir. uj 2 — (u/loq) 2 is now 
eigenvalue of $ = (& 5 /16)<I>. For practical purposes, the 
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FIG. 1: a) First Brillouin zone for a triangular lattice. The 
corresponding irreducible zone is delimited by the points F, 
J and X. b) Dispersion curves for phonon modes in the large 
r s limit (crystal). Each branch corresponds to a particular 
polarization. The wavevector q is taken on the boundary of 
the irreducible Brillouin zone along the path V — > J — > X — > 

r. 



kernel ^(q) is not computed by direct summation but we 
use a variant of the Ewald summation technique instead, 
see appendix 

Correction to the ground state energy with respect to 
the classical case (r s — > +00) stems perturbatively from 
the zero-point energies of these phononic modes. Explic- 
itly, the ground state energy per particle is given by 



R, 



E 

i=i,2 



uz v bz 2 



where j = 1,2 stands for the two phonon branches and 
V B z = 8ir 2 /V3b 2 is the area of the Brillouin zone (BZ). 
In reduced units it takes the form 



E 




d 2 q 



IBZ 



where the summation is performed over the irreducible 
Brillouin zone (IBZ), see Fig. [T^,)). The numerical eval- 
uation gives 



E = 1.59702 



2.02987 



O 



(13) 



At a not too large value of r s — 50 and for N = 50 
particles the GFMC simulation gives Eq FMC = 1-893 
while the aforementioned expansion provides E = 1.884, 
i.e. within 0.5% of the exact result. 



The ground state wavefunction is a simple product of 
all Gaussian ground states of the phonon modes. It can 
be written as 



$g({u,}) =Me 



5Ti = l,2 f B 



with A/" a normalization factor. Vectors u q are Fourier 
transforms of the original variables {u^} and eigenvectors 
of the system ©. The reduced onebody wavefunction 
is obtained by integrating all atomic positions but one. 
This leads to the following anisotropic Gaussian wave- 
function 

*G,i(aj, v) = M^-^^+^y^/ 2 , 

with typical sizes Xq and respectively in the x and 
y directions. Let (a(q),/3(q)) and (— /3(q), a(q)) be the 
normalized (a(q) 2 + /3(q) 2 = 1) eigenvectors of Eq. 
we find 



2 _ 3^3,9 



6 9/2 



d 2 q /a(q) 2 , P(qf 



IBZ 



(2tt) 2 V^i(q) Wa(q) 



and the same expression holds for yo where a(q) and /3(q) 
are interchanged. A numerical evaluation yields xq 



1.3613/rJ /4 and y Q = 0.9509/r^ /4 in units of 



1/4 



IV. FIRST ORDER NATURE OF THE 
QUANTUM MELTING 

In this section, we use the quantum Monte-Carlo cal- 
culations in order to determine the order of the transi- 
tion. We find strong numerical evidences of a first order 
transition. 

Close to a first order transition, the GFMC energies are 
expected to depend on the GWF \^g)- To illustrate this 
point qualitatively, let us expand \^g) on the eigenbasis 
of the Hamiltonian, 



|*Gf(a)> oc |*bec) +c(a)|* 



cry/ 



(14) 



where I^bec) (l^cry)) stands for the BEC (crystal) 
eigenstate and c(a) measures the relative projection of 
the GWF on the two competing eigenstates (a more care- 
ful treatment would also include the low lying excitations 
of the two phases). After projection, one gets a state 
I*) = e-^|* G (a)), 



I*) oc |*bec> +c{a)e 



-/3(Bbe 



(15) 



where -Ebec (-^cry) stands for the BEC (crystal) energy. 
For very large imaginary time /3, the exponential is very 
large (or small), and 1^) converges to whichever of the 
BEC or crystal state is more stable. Real simulations 
however are done with finite values of (3 (typically 30/i 
in our case). As one approaches the transition, -Ebec — 
E cj: y vanishes, hence the exponential term tends to one, 
and one is left with a (incoherent) superposition |\&) oc 
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FIG. 2: GFMC calculation of the energy per particle for var- 
ious values of the symmetry breaking parameter a = (cir- 
cles), a = 0.2 (empty squares), a — 0.65 (diamonds) and 
a — 0.75 (empty diamonds). The transition is emphasized by 
plotting Bgfmc - 1.59702 - 2.02987/^/rI - 0.85/r s as a func- 
tion of r s (the last number here, 0.85, is a fitted parameter). 
Above the critical r s ~ 27± 1, the symmetry breaking state is 
energetically favored. This plots also demonstrates that the 
transition is first order. The system contains N = 72 bosons 
in 48 x 84 sites. 



I^bec) +c(a)|vE'cry) whose properties are sensitive to the 
choice of the variational parameter a of the GWF. 

In Fig[2j we present the GFMC energies as a func- 
tion of r s for various choices of the symmetry break- 
ing parameter a. Those data support the previous pic- 
ture: concentrating at first on the a = ("BEC" like 
state, with c(0) <C 1) and a = 0.75 ("crystal" like state, 
with c(0.75) 3> 1) curves, we find a crossing around 
r s as 27 above which the crystal becomes more stable 
than the BEC. When the symmetry is only weakly bro- 
ken (a = 0.2), the GFMC result is closed to the BEC 
state close to the transition. As one increases r s fur- 
ther, the energy difference £bec — E CYy increases and 
the exponential in Eq. (fT5|) eventually compensates the 
smallness of c(0.2) and |^) converges toward the actual 
ground state |^ C ry)- 

To further characterize the transition, we now look at 
the crystalline order parameter. Introducing the static 
density-density correlation function g{r) (that roughly 
measures the probability of finding a particle on point r 
knowing that there is one at point 0), 



L X Ly 



E 



N(N + 1) 

h 



r+h h h r+h/ 



(16) 



the order parameter is usually defined as the Fourier 
transform of g(r) taken at a k vector belonging to the 
reciprocal lattice of the expected crystal. A typical ex- 
ample of g(r) in the BEC (crystal) case can be found in 
Fig. 2k (Fig. Hb). Here we define the crystalline order 
parameter G(r s ) (in an equivalent but more convenient 
way) as the average value of g(r) at the peaks corre- 




FIG. 3: Crystal order parameter G(r 3 ) for three value of the 
symmetry breaking parameter a = (circles), a — 0.4 (empty 
squares) and a = 0.8 (diamonds). The system contains N = 
72 bosons in 48 x 84 sites and G(r s ) has been evaluated at 
/3 = 30. The vertical dashed line indicate the region where the 
results are very sensitive to the choice of the wave-function. 



sponding to the classical position of the crystal (minus 
its average value). In FigJ3J we plot G(r s ) for non broken 
(a = 0), broken (a = 0.8) and partially broken (a = 0.4) 
symmetry. The non broken and broken symmetry case 
have respectively zero and non zero G, hence G is discon- 
tinuous at the transition which should therefore be first 
order. The case with partially broken symmetry shows 
the expected crossover near the transition point. 

Far from the transition point, the results are supposed 
to be independent of the GWF. We checked that this is 
indeed the case in our calculation: Fig. 0] presents the re- 
sults before (Variational Monte-Carlo VMC, on the left) 
and after (RQMC, on the right) projection on the ground 
state. One can see that although the starting point has 
no broken symmetry (a = 0), the RQMC (or GFMC) 
algorithm is able to converge to the crystal ground state 
nevertheless. This result was obtained for a rather small 
system of N = 32 particles however, and for a larger sys- 
tem of N = 72 bosons, a small amount of a was needed 
in order to get the correct ground state. 

The above results are consistent with those obtained in 
Refi 19 ' 20 . In Ref. [l^, the authors use a finite temperature 
algorithm and find the transition at r s — 32 ± 7. Indeed, 
in our unit, they work at a temperature of the order of 
5. 10~ 3 , i.e. several time the difference of energy (per 
particle) between the BEC and the crystal phase near 
the transition, see fig. |U These explains the relatively 
important error bar on their critical value of r s and the 
fact that their superfluid fraction switches back and forth 
between and 1 close to the transition. In Refi^i, the 
authors use the DMC technique, similar to ours. The 
GWF that they are using to describe the crystal is not 
symmetric with respect to the interchange of the bosons, 
and hence describe a crystal of distinguishable particles. 
However they find a critical r s = 30 which is consistent 
with ours (we found upon decreasing the filling factor v 
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FIG. 4: Density-density correlation function g(f) at r s = 50, 
well above the melting point. The calculations are done for 32 
particles in a 32 x 56 grid with the trial wave-function for the 
condensate, a = 0. On the left, the VMC results (a) show no 
sign of broken symmetry while after projection on the ground 
state, the RQMC results (b) for f3 — 20 are in a crystal state. 
The gray scale goes from zero (white) to 4 (black). 



a small increase of the critical value of r s ). 

We now proceed with the Maxwell construction, in or- 
der to find the region where phase separation is to be 
expected. In the rest of this section, we measure ener- 
gies in unit of Eg = Cdd/ipa^) and not Eg as in the 
rest of the paper since Eb depends on r s . Denning the 



density n 



such a phase separation can occur be- 



tween the two values of the density n\ and n 2 satisfying 
[E(ni) - £(na)]/(fti - n 2 ) = dE/dn{n{) = dE/dn(n 2 ). 
Close to the transition, we find that a good parameteri- 
zation of the energy reads, 



E w a n 3/2 +a 1/2 n 5/ ' i +a 1 n+r]6(Vn-r* s 

J (17) 

where a ~ 1.592, a 1/2 w 2.03, a x w 0.85, r) w 5 10~ 4 , 
r* w 27, and 9{x) is the Heaviside function discriminat- 
ing the two phases. Note that the coefficient r\ can be 
directly extracted from Fig. [5J We find that the region 
where phase separation can occur is very small, given by 
a width Sr s = a(r*) 2 /(3ao) around the transition point 
r*. For our parameters, it translates into Sr s « 0.04 
which is very small. 



= 5/4. 



-3/2 



V. ON THE PRESENCE OF A SUPERSOLID 
PHASE 

Once the presence of a first order transition has been 
established, a natural question is whether the superfluid 
fraction vanishes right at the transition or if there is a 
region where both the superfluid fraction and the order 
parameter G are non zero simultaneously. Such a region 



would be a supersolid phase and has been searched for 
for many years. However, to our knowledge, the only 
model for which it has actually been found are bosons on 
a lattic e) 27 ' 28 ! 29 at large filling factors (i.e. systems similar 
to our but in the limit v 1 not v <C 1) where the crystal 
is somehow stabilized by the underlying grid. When the 
crystal is uniquely due to the repulsive interaction (as in 
Helium 4 or in the system considered here), no convincing 
evidence could be gathered so far. 

There are two alternative way to probe the superfluid 
character of the system. The direct way consists in cal- 
culating directly the superfluid fraction ps of the system. 
It is denned as the fraction of the bosons that does not 
follow when the frame is put (adiabatically) into motion 
(rotation). It is the bosonic equivalent of the curvature 
of the energy with respect to an Aharonov flux for elec- 
trons. For quantum Monte-Carlo calculations, it is can 
be conveniently calculated using the diffusive motion of 
the center of mass of the system (in imaginary time) as 



PS 



myi- 



tis) 



where (3 is the imaginary time and R 2 {(3) is the second 
moment of the center of mass in one direction (say x). 
Such a formula can be easily implemented, but the re- 
sults should be taken with care. For r s < 27 (BEC state) 
we find that ps ~ 1- For large r s , we find that ps ~ 0. 
However, for r s > 27 but close to the transition, we could 
not get a reliable measure of ps and hence cannot con- 
clude about the presence of a supersolid in this model. 
This is due to a combination of two difficulties: (i) close 
to the transition, "metastable state" can enter into the 
GWF, hence artificially enhancing ps (depending on the 
choice of (3, see the discussion in section ll*V|) . (ii) In the 
crystal phase, the superfluid fraction is given by rather 
rare events that contain "dynamical defects" 7 . To il- 
lustrate the difficulty of the calculation, let us look at 
a small system of N — 8 particles, where using brute 
force simulations, reliable results could be obtained for 
any choice of a. Using a very large number of walkers 
(10 6 ), and rather long simulations (/3 = 80/i), we find 
in Fig|5]a superfluid fraction ps ~ 0.82 for all values of 
a. However, for a = 0.5 the correct superfluid fraction is 
only recovered after a (large) simulation time (j3 > 20). 
A careful analysis shows that for (3 < 20, the energy 
is higher than the ground state energy by a very small 
amount (of the order of 10~ 4 ) usually undetectable in 
reasonable calculations. In the inset of Fig[5] (a zoom 
of the main figure which would correspond to a typical 
run for a larger system) we find that the datas are per- 
fectly fitted by R 2 {(3)/(2N) = 0.41/3+1. 24(1 -e"' 3 / 3 - 4 ) (a 
usual form) leading to the much smaller (and incorrect) 
ps ~ 0.41. Hence, we find that it is very important for 
practical calculations to check both the convergence of 
the energy and the dependence of ps on the GWF very 
carefully. A poor choice of the GWF will lead to a slow 
convergence of the energy (see later in this section) and 
to incorrect measure of the superfluid density. 




FIG. 5: (color online). Diffusion of the center of mass R 2 ([3) 
for N = 8 particles in 16 x 28 sites at r s — 200 using 10 6 walk- 
ers. The various lines correspond to a — (straight), a = 0.2 
(dotted) and a — 0.5 (dashed). The thin red (blue dashed) 
lines are Et y = 0.82/3 ( y = 0.41/3 + 1.24(1 - e~ p/3A ) ) corre- 
sponding to a superfluid fraction of 0.82 and 0.41 respectively. 
Inset: zoom of the main figure. 



A second way to find a supersolid is through non diag- 
onal long range order in the one body density matrix 



9i(r,r) = — ^ 



etc, 



By definition, its largest eigenvalue gives the conden- 
sate fraction no- For homogeneous systems it leads to 
n = (l/L x L y ) J2 r 5i( r + h, h)- Alternatively, n is also 
given by the long-range asymptotic |r — r'| — > +oo of 
gi(r',r) and the two definitions coincide in the thermo- 
dynamical limit. In our case, we reduce significantly fi- 
nite size effects by excluding the short range part of g\ 
with an arbitrary length ro, 



n 



— J2 5i(h + r,h), 

Vr 



r|>r 



where V ro is the system volume excluding a disk of radius 
r . In contrast to density-density correlations, GFMC 
only provides access to mixed estimators for the con- 
densate fraction since gi(r',r) is a non-local quantity. 
Results then depend on the quality of the guiding wave- 
function compared to the exact wavefunction. For small 
r s , interactions are weak and our variational Bijl-Jastrow 
form ^ is quite a good approximation. In that case, 
the bias due to mixed estimators can be further reduced 
by computing 2n vIX — n /MC where n /MC is the varia- 
tional Monte Carlo result using the same guiding wave- 
function. The situation is rather different close to the 
quantum melting point since the Bijl-Jastrow does not 
describe the strong correlations between particles. For 
this case, the mixed estimator for g\ leads to results that 
can not be trusted. To illustrate this point, we plot the 
condensate fraction for r s = 27 and N = 50 particles 
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FIG. 6: Condensate fraction no as a function of the variational 
parameter a, see Eq. {5b. The three lines correspond to a 
VMC calculation with tiq MC , a GFMC calculation which gives 
a mixed estimator no* and an extrapolated result 2n VIX — 
n yMC This i as t quantity is reduced by a factor 3 between 
a = 0.5 and a = 0.8. 



and for different values of a in Fig. O The strong de- 



pendence on a even for 2n, 



MX 




,VMC 
l 



clearly rules out 



the relevance of mixed estimators to compute condensate 
fractions. Note that the results presented Fig.[S]are quan- 
titatively only weakly modified after a finite size scaling 
analysis. In particular, it is not possible to conclude from 
these calculations on the presence of a supersolid with a 
non vanishing condensate fraction. 

The conclusion of this part is that, although it is not 
difficult to calculate the energy of the system, a better 
GWF is actually needed to get a reliable calculation of 
superfluid properties. Near the transition point, our un- 
derstanding of the physics, and hence the corresponding 
GWF, is not sufficient. The quality of the GWF is usu- 
ally measured using either its energy (difference to the 
ground state) or the variance of its energy (that vanishes 
for eigenstates). For practical calculations however, some 
GWF can have a decent variational energy, but converge 
very slowly toward the ground state, hence being poor 
choice in practice. Following the recent proposal of two 
of us^, we measure the quality of the GWF, with its 
overlap f2 with the actual ground state of the system. 
Defining 



_ -kN _ 



(*gI*o 



(*o|*o)(*g|*c 



(19) 



the parameter k is given by22, 

dp [E{fi) - E{p = +oo)] (20) 



where E(f3) is the energy (as a function of imaginary time 
and per particle) of the system while E(f3 = +oo) is the 
ground state energy (per particle). Hence k measures 
the area below the curve E{0) and is in this sense a good 
measure of the practical use of the GWF (a slow con- 
verging GFW will have a large k for instance). A study 
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FIG. 7: Parameter k measuring the overlap between the guid- 
ing function and the actual ground state of the system. k(t s ) 
at a — (squares) and a — 0.5 (circles). The datas are taken 
for 32 bosons in a 32 x 56 grid. 



of k shown in Fig. [7] reveals two things: (i) The broken 
symmetry GWF (a — 0.5) actually gets better than the 
BEC (a — 0) around the transition point r s = 27. (ii) 
However, in the broken symmetry phase, n still increases 
with r s , indicating that we describe poorly the crystal 
close to the transition. 



VI. DISCUSSION 

We conclude this paper by discussing how the above 
results would apply to actual experiments using ultra- 
cold dipolar molecular systems. The system we have in 
mind could be for instance, RbCs molecules confined by a 
two-dimensional trap and where electric dipole moments 
are induced by a perpendicular electric field. Those 
types of systems typically have three different experi- 
mental limitations that constrain the range of parame- 
ters (r s ) than could actually be observed in practice, (i) 
First, the dipolar forces that can be induced is limited 
to Cdd = 7(eao) 2 /(47re ) (eo dielectric constant, a Bohr 
radius) with 7 < 7 cxp (typically 7 cxp ps 0.25 for RbCsi^). 
(ii) The maximum density is also limited and £ > £ cxp 
with a typical value of Z exp ps 60nm. It means that r s 
is limited by r s < (£* B /£ C x P )- Noting £* B = j£ , we get 
r s < 7(4/4x P ) where £ = M (ea ) 2 / (ft 2 47re ) ps 20/xm 
for RbCs. (iii) the temperature is limited to T exp of 
the order of T exp ~ lOOnK. In order to observe the 
physics discussed here T exp should be a fraction e (say 
e fa 1%) of E B . This imposes, r s > (T cxp /eT ) 1/3 7 2/3 
with fcT = (ea ) 2 /(8ne e%). For RbCs, T w 1.45nK. It 
is an experimental challenge to actually produce ultra- 
cold heteronuclear molecules in their vibrational and ro- 
tational ground state. However we stress here that, once 
these molecules are produced, the quantum melting tran- 



FIG. 8: Experimental accessible range of parameters in the 
(7, r s ) plane for a 2D system of RbCs dipolar molecules. The 
solid lines indicate the three experimental constraints dis- 
cussed in the text, while the dashed line stands for the critical 
value of r s m 27. The white region should be accessible to ex- 
periments. 

sition itself should be within reach. Putting these three 
above conditions together, we find the allowed window of 
parameters shown in Fig. [8] 
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APPENDIX A: EWALD SUMMATION 

Finite size effects can be reduced by adding periodic 
images to the original simulation box forming a lattice 
pattern with rectangular boxes of size L x x L y . Each 
atom in the original box has images in all other boxes. 
The simple 1 /r 3 dipole interaction between atoms is thus 
replaced by an effective interaction taking all images into 
account, or 

^-r'> = SI |r _ r ,; RMr (AD 

where Rm = (nL x ,mL y ) denote all vectors of a rectan- 
gular lattice. Since the summation has a poor conver- 
gence, we extend the well-known Ewald summation trick 
to this problem. The y-integral in the identity 

^ = A/ d yy 2 e-^ 2 , (A2) 
|r| 3 Jo 

is splitted in two parts J^°° = + / + °° and V(r) is 
splitted correspondingly V(r) = V < (r) + V > (r). yo is 



9 



a parameter that can be chosen arbitrarily. The second 
part for y > y converges exponentially fast when the 
summation over positions is performed. One finds 

V > {r)=y$ J2^(yo\r + R M \). 

We have introduced the set of functions 
2 r+oo 



<p n (x) = -== / dtt n e~ tx 



(A3) 



which can be expressed as rational functions of Erf (a;), 
e~ x and powers of x. For the case y < yo, we define the 
function F(r,y) = Y^b. e V ' t+Rm ' periodic over the 
rectangular lattice for fixed y. Thus we write F(r, y) as a 
summation over its Fourier components on the reciprocal 
lattice 



F(r,y)=J2e iK » 



Km 



7r e Ml 



L x L y y 2 



with reciprocal lattice vectors Km = (2im/L x ,2Trm/L y ), 
and the relation 

4 f yo 

V<(r) = n = dyy 2 F(r,y), 

V 7 *" JO 

yields again an exponentially fast converging expression 
for V^ < (r). After some straightforward calculations and 
using definition (|A3|) for functions <*p n (x) we finally arrive 
at 



V(r) =Vb i J2 Vi/a(W) l r + R : 

Rm 

1/0 V"^ , v s f\ K M\ 

+ j—j- C0S ( K A/ ' rV_3/2 



2yo 



(A4) 



the effective interaction between atoms in the simulation 
box. Moreover, subtracting l/|r| 3 from the right-hand- 



side of Eq. (| A4[) and letting r — > leads to the constant 
A in Eq. ©, 



A = — 



R 



Ml 



[Km] 
2y 



which accounts for the self-interaction between each atom 
and its own images. 

The same procedure can be applied to calculate effi- 
ciently the kernel (flQ|) describing phonons for the crystal 
phase. Again we start with Eq. (|A2j) . split the y- integral 
and $ = <! < +<i> > accordingly. The second part converges 
very rapidly when the summation over the positions Rj 
of the triangular lattice is performed, one finds 



■<P5/2 



\TLi 



, l R ^l 



(l-cos(q- Rj)), 



with functions (|A3|) . For convenience, we choose yo = 1/b 



where b — y 2/V3 is the triangular lattice spacing. The 
first part is computed more efficiently by going to Fourier 
space in a similar way as above. The final result reads 



'b\G s +q\\ ,_ ,_ fb\Gj\ 



where the summation is carried out over all vectors Gj 
on the hexagonal reciprocal lattice with an exponentially 
fast convergence. 
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